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ABSTRACT 

In this thesis single-degree-of-freedom torsional airfoil flutter is investigated 
using an incompressible potential flow code, a compressible inviscid Euler code and a 
compressible viscous Navier-Stokes code. It is found that the classical linearized 
incompressible and compressible flow theories yield unconservative flutter estimates. 
The computations based on the non-linear codes show for NACA 0006, NACA 0009, 
NACA 0012 and NACA 0015 airfoils that the regions of torsional flutter instability 
increase as the airfoil thickness and the flight Mach number is increased. On the other 
hand, the comparison of the flutter boundaries computed with the viscous Navier-Stokes 
code versus the inviscid Euler code shows that the effect of viscosity is stabilizing. Also, 
the computed flutter boundaries display the effect of pitch axis location on flutter. Axis 
locations in the range between half a chord upstream of the leading edge of the airfoil and 
the leading edge are most prone to induce flutter. Axis locations downstream of the 
quarter chord are flutter free. 


V 





TABLE OF CONTENTS 


l. INTRODUCTION. 1 

A. BACKGROUND. l 

B. OBJECTIVE.2 

n. GOVERNING EQUATIONS.3 

A. CONTINUITY EQUATION.3 

B. MOMENTUM EQUATION.•..3 

C. ENERGY EQUATION.4 

D. VECTOR FORM OF FLOW EQUATIONS.5 

E. TURBULENCE MODEL.6 

m. NUMERICAL FLOW SOLUTIONS. 9 

A. INCOMPRESSIBLE INVISCID FLOW.,.. 9, 

B. COMPRESSIBLE FLOW.10 

1. Grid Generation by Algebraic Mapping...12 

2. Numerical Implementation of Algebraic Mapping.13 

3. CFD Techniques.14 

a. Discretisation process example.15 

b. Convergence..17 

4. Navier - Stokes Solution. 17 

IV. AEROELASTIC ANALYSIS.‘.21 

A. FLUTTER - GENERAL.21 

B. THE PHENOMENON OF FLUTTER.22 

C. NONDIMENSIONAL PARAMETERS.24 

D. STIFFNESS CRITERIA.27 

E. VERTICAL TRANSLATORY OSCILLATION.27 

F. ONE DEGREE OF FREEDOM FLUTTER (TORSIONAL).30 

vii 




























1. Physical Explanation.30 

2. Flutter Analysis.33 

a. Time domain analysis.34 

b. Frequency domain analysis.36 

G. TWO DEGREES OF FREEDOM FLUTTER.37 

1. Physical Explanation.37 

2. Equations of motion.40 

a. Time domain analysis.40 

b. Frequency domain analysis.42 

H. NS - UPOT REDUCED FREQUENCY COMPATBILITY.43 

I. FLUTTER PREDICTION WITH UPOT CODE. 44 

V. RESULTS.47 

A. GENERAL 47 

B. GRID GENERATION. 47 

C. INPUT DATA DESCRIPTION . 49 

D. CALCULATION PROCEDURE... 52 

E. NON - VISCOUS (EULER) CALCULATION.53 

F. EFFECT OF PIVOT POINT ON TORSIONAL FLUTTER.63 

G. EFFECT OF MACH NUMBER ON TORSIONAL FLUTTER.64 

H. EFFECT OF AIRFOIL THICKNESS.67 

I. VISCOUS (NAVIER-STOKES) CASE CALCULATION...70 

J. EULER-NS RESULTS COMPARISON.76 

K. PANEL - EULER CODE COMPARISON.83 

VI. SUMARY. 85 

Vn. RECOMMENDATIONS... 87 

LIST OF REFERENCES .89 


viii 


INITIAL DISTRIBUTION LIST 


91 





























LIST OF FIGURES 


Figure 3-1 Incompressible Flow wake model.10 

Figure 3-2 Overview of Computational Fluid Dynamics.11 

Figure 3-3 Numerical Implementation of Algebraic Mapping.14 

Figure 3-4 Overview of the computational solution procedure.15 

Figure 4-1 Vector diagram of lift in vertical translation oscillation. 29 

Figure 4-2 Oscillating airfoil about a pivot point inf front of Vi chord.31 

Figure 4-3 Incompressible flow wake visualization (from UPOT).32 

Figure 4-4 Rigid, symmetrical airfoil restrained to rotate about L.E.33 

Figure 4-5 Variation with reduced frequency k of the real and imaginary parts of the dimensionless 
aerodynamic moment m^ due to pitching of an airfoil about its leading edge in 

incompressible flow.37 

Figure 4-6 Two degrees of freedom airfoil motion - phase difference 90®.38 

Figure 4-7 Torsion - Bending diagrams- phase difference 90®.38 

Figure 4-8 Two degrees of freedom airfoil motion - phase difference 0®.39 

Figure 4-9 Torsion - Bending diagrams- phase difference 180®.39 

Figure 4-10 Aerodynamic moment coefficient graph for low k values.. 44 

Figure 4-11 Aerodynamic moment coefficient graph for high k values.45 

Figure 4-12 Aerodynamic moment coefficient graph for the critical flutter k value.45 

Figure 5-1 C-grid for NACA 0015 (non-viscous flow) 48 

Figure 5-2 NACA 0015 C-grid details (non-viscous flow). 49 

Figure 5-3 Steady state solution for NACA 0015 M=0.7 (Euler case).54 

Figure 5-4 Unsteady case solution - Ka=0.36 (Euler case).55 

Figure 5-5 Time rate of change of AOA amplitude for Ka;=0.36 (Euler case).56 

Figure 5-6 Unsteady case solution - Ka=0.25 (Euler case).56 

Figure 5-7 Time rate of change of AOA amplitude for Ka =0.25 (Euler case). 57 

Figure 5-8 Curvefit of Ka vs time rate of change of AOA amplitude (Euler case).58 

Figure 5-9 Unsteady pase solution - Ka=0.324 Critical flutter case (Euler)...58 

Figure 5-10 Time rate of change of AOA amplitude for Ka =0.324 Critical flutter (Euler case).59 

Figure 5-11 Pressure Contours around NACA 0015 for a whole cycle of oscillation M=0.7 Xp=0.0.61 

Figure 5-12 Effect of Pivot Point on reduced frequency for NACA 0015 and M=0.7 (Euler case).63 

Figure 5-13 Effect of Pivot Point on critical k^ for NACA 0015 and M=0.7, ia=100 (Euler case).64 

Figure 5-14 Effect of Mach Number and Pivot Point for NACA 0015 (Euler case).65 

ix 































Figure 5-15 Reduced frequency variation with Mach number for Pivot Points -0.85<Xp<0.15 

(NACA 0015).65 

Figure 5-16 Effect of Mach number and pitch axis location on torsional flutter of NACA 0015 airfoil, 

ia=100...66 

Figure 5-17 Variation of flutter speed with Mach Number for pivot points -0.85<Xp<0.15, Ia=100.67 

Figure 5-18 Effect of airfoil thickness at M=0.7.68 

Figure 5-19 Effect of airfoil thickness at M=0.7.68 

Figure 5-20 Effect of airfoil thickness on torsional flutter at M=0.7, ia=100.69 

Figure 5-21 Yariation of flutter speed with airfoil thickness for various pivot points -0.85<Xp<0.15, 

icFlOO. 70 

Figure 5-22 Steady state soln for NACA 0015 M=0.7 (N-S case).71 

Figure 5-23 Unsteady case solution - Ka=0.34 (N-S case)..72 

. Figure 5-24 Time rate of change of AOA amplitude for Ka =0.34 (N-S case).73 

Figure 5-25 Unsteady case solution - Ka=0.20 (N-S case).73 

Figure 5-26 Time rate of change of AOA amplitude for Ka =0.20 (N-S case).74 

Figure 5-27 Curvefit of Ka vs time rate of change of AOA amplitude (N-S case).75 

Figure 5-28 Unsteady case solution - Ka=0.252 Critical flutter case (N-S).75 

Figure 5-29 Time rate of change of AOA amplitude for Ka=324. Critical flutter (N-S case).76 

Figure 5-30 Euler - NS reduced frequency results for NACA 0015 and M=0,7.77 

Figitfe 5-31 Euler - NS reduced natural pitching frequency results for NACA 0015 and M=0.7, 

ia=100...^.78 

Figure 5-32 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (Euler).79 

Figure 5-33 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (Euler-detail).. 79 

Figure 5-34 Mach contours for NACA 0015 M=0.7 steady AOA (Euler).80 

Figure 5-35 Mach contours for NACA 0015 M=0.7 steady AOA (Euler-detail).80 

Figure 5-36 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (N-S).81 

Figure 5-37 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (N-S-detail).81 

Figure 5-38 Mach contours for NACA 0015 M=0.7 steady AOA (N-S).82 

Figure 5-39 Mach contours for NACA 0015 M=0.7 steady AOA (N-S-detail).82 

Figure 5-40 Low subsonic Euler code results comparison with UPOT results.83 

Figure 5-41 % differences between Euler and UPOT results (w.r.t M=0.1 results).84 


X 






























LIST OF TABLES 

Table 5-1 Steady state solution for NACA 0015 M=0.7 (Euler case).53 

Table 5-2 Regression equation slope values for various ka values.57 

Table 5-3 Average period T calculation for NACA 0015 flutter - Ka=0.324.60 

Table 5-4 Percentage Increase of k with Mach number and Pivot Points.66 

Table 5-5 Percentage increase of k with airfoil thickness and pivot point location.69 

Table 5-6 Steady state.solution for NACA 0015 M=0.7 (N-S case). 71 

Table 5-7 Regression equation slope values for various kot values. 74 

Table 5-8 Numerical differences of reduced frequency between Euler and N-S calculations for 

NACA 0015 and M=0.7.;.77 


xi 













LIST OF SYMBOLS 


b 

c 

Cl 

Cm 

Cp 

e 

E 

f 

h 

la 

ia 

J 

k 

ka 

h 

Ka 

Kh 

L 

m 

Ma 

P 

q 

Re 


free-stream speed of sound 

airfoil half chord 

airfoil chord length 

lift coefficient per unit span 

pitching moment coefficient per unit span 

pressure coefficient 

total energy per unit volume 

energy 

frequency 

vertical displacement in plunge defined positive downward 
dimensional moment of inertia 

non-dimensional moment of inertia i„ = —— 

TtpyiA 

transformation Jacobian flux of momentum 
reduced frequency, iTtfctU^ 
reduced natural pitching frequency 
reduced natural plunging frequency 
spring constant for pitching 
spring constant for plunging 
lift per unit span 
mass 

free-stream' Mach number 
aerodynamic moment per unit span 
pressure 

flow field vector 
Reynolds number 

xiii 




t 

T 

u,w 

Vf 

W 

Xq 

Xp 

x,z 

a 

r 

A 

p 

p 

T 

0 

0)a 

COh 


Static moment (= x^m) 
time 

period of oscillation 
free-stream velocity 

velocity components at x,z direction respectively 

reduced flutter velocity, = 1 / 

work 

center of mass to elastic axis location distance 

leading edge to elastic axis location distance 

coordinates in physical plane 

airfoil angle of attack 

circulation 

vorticity strength 

wavelength of oscillation 

viscosity 

coordinates in transformed plane 
density 

non-dimensionalized time 
velocity potential 

Uncoupled natural torsional frequency, 
uncoupled natural pitching frequency, ^K^Jm 


4 


Operators 


0 

o' 


differentiation with respect to time 
differentiation with respect to non-dimensional time 


XIV 



A 


incremental change 
del operator 


V 


XV 





ACKNOWLEDGEMENT 


I would like to express my sincere appreciation to my thesis advisor. Professor 
Platzer. His experience and deep knowledge in aeroelasicity proved to be extremely 
beneficial for me. He has his own way to make the student understand the essence of this 
science. The study of aeroelasticity phenomena is an extremely difficult job but becomes 
fun when working with Professor Platzer. I owe most of this work to him. 

Many thanks to Dr. K.Jones. None of this work could have been completed 
without his help on the computer codes and the use of software. He spent a lot of his 
valuable time in order to show me how I should proceed. He was always there whenever 
I needed him. 

Finally I would like to thank Dr. S. Weber. He was the one who helped me a lot in 
understanding the way the codes work. He was also the one who advised me what to do 
when everything seenred to be stuck and the “no way out” seemed to be unavoidable. 


XVll 



I. INTRODUCTJON 


A. BACKGROUND 

The phenomenon of ahfoil flutter was first encountered on World War I airplanes. 
The danger to flight safety posed by this phenomenon stimulated the development of 
quite sophisticated flutter analysis and testing methods in the ensuing years. References 
1,3,4,6 and 7 give good reviews of these developments, where it is also explained that the 
flutter phenomenon is caused by the interaction of the aerodynamic, inertia and elastic 
forces. Until very recently, important simplifying assumptions had to be introduced in 
order to make the problem mathematically tractable. Linearized aerodynamic analysis 
methods were among the most important simplifications used in the past. For the analysis 
of low speed flutter phenomena, Theodorsen developed an incompressible flow theory 
for oscillating flat plates, which was then extended to compressible subsonic flow. This 
incompressible oscillatory flat plate theory remains the standard aerodynamic analysis 
tool for flutter calculations. A second important simplification was also introduced by 
Theodorsen when he proposed to perform the flutter analysis in the frequency domain. In 
this approach, the aerod 5 mamic forces need to be computed only for the special case 
where the airfoil is assumed to execute a purely harmonic oscillation at constant small 
amplitude. 

Flutter analyses based on linearized aerodynamic theory in the frequency domain 
have serious shortcomings. Thin airfoil (flat plate) theory makes it impossible to account 
for the effect of airfoil geometry on flutter. Frequency-domain methods provide little 
information about the physics of the flutter problem because the decay of the airfoil 
motion below the critical flutter speed and the divergence of the motion above this speed 
cannot be obtained as part of frequency domain flutter analysis. 

The recent rapid advances in computational fluid dynamics (CFD), however, 
made it possible to replace linearized aerodynamic methods by solutions based on the full 
non-linear Euler or Navier-Stokes equations for inviscid or viscous compressible flow. 
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As an additional advantage, these modem CFD methods integrate the governing flow 
equations in time. Therefore, they can easily and naturally be combined with the 
equations of motion of the airfoil. This leads to the following specific objective for the 
work attempted in this investigation. 


B. OBJECTIVE 

A two-dimensional compressible Euler and Navier-Stokes flow solver is coupled 
with a one degree-of-freedom stmctural model for the time domain analysis of an airfoil 
which is free to oscillate about a specific pivot point. The effect of airfoil thickness on 
torsional flutter is to be determined as a function of pivot location and subsonic Mach 
number. Also, viscous flow effects on flutter are to be determined by comparing Euler 
and Navier-Stokes computed flutter boundaries. Furthermore, at low Mach numbers the 
Euler computed flutter boundaries are to be compared with the flutter boundaries 
computed with the incompressible unsteady panel code UPOT. 
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II. GOVERNING EQUATIONS 


The equations that describe the compressible viscous fluid flow around a body are 
the continuity, the momentum and the energy equations. The flow around the body can be 
computed by the simultaneous solution of these equations. The conservation-law and the 
vector form of the compressible Reynolds-averaged Navier-Stokes equation is presented 
in this chapter. A detailed derivation can be found in [Ref 8]. 


A. CONTINUITY EQUATION 


The continuity equation expresses the conservation-of-mass law applied to a 
fluid passing through a control volume fixed in space 

^ + V(py) = 0 (2.1) 

at 

where p is the fluid density and V is the fluid velocity. Eq.(2.1) states that the net mass 
flux through a control volume bounding surface must be equal to the time rate of change 
of the mass inside the control volume. For two-dimensional Cartesian flow this equation 
reads 

+ ( 2 . 2 ) 

where u and w are velocity components along the x and z directions, respectively. 


B. MOMENTUM EQUATION 

The momentum equation expresses Newton’s second law as applied to a fluid 
element flowing relative to a space-fixed coordinate system. The x and z direction 
momentum equations are: 
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(2.3) 


^^^+V(/7yy)=// + vn, 

at 

The first term in Eq.(2.3) represents the time rate of change of momentum per unit 
volume in the control volume. The second term represents the momentum net flux 
through the bounding surface of the control volume. Also, f is the body force per unit 

volume and y the stress tensor given by: 


dui 2^ 
dxj 3 dx 


k J 


(2.4) 


where i,j,k = 1,2,3 and 6ij is the Kronecker delta. 

By substituting Eq.(2.4) into (2.3) for flow in a two-dimensional Cartesian 
coordinate system we obtain: 
dpu dpu dpu . dp d 



d 


^ dw ^ du \ 

2/3//(2— 

dx dz _ 

^ dz 


^dx dz)_ 


(2.5) 


dpw dpw dpw . dp d 


2/3// 


( 2 ^^ ^ ^ 

di dx. 




dx 


dw • du 


These equations are known as the Navier-Stokes equations. 


C. ENERGY EQUATION 

The energy equation is derived by applying the first law of thermodynamics (rate 
of change of energy = net heat flux into particle + rate of work done on particle). 

p(^f u-\r f w) + ^{eu + pu-uT^ -WT^)+ (2.6) 

3/ dt ax 

—(ew + pw-wt^ -uT^)-0 
■ dz 

where e is the total energy per unit volume. 
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D. VECTOR FORM OF FLOW EQUATIONS 

The above equations can be rewritten in non-dimensionalized vector form as 


where 


Q = 


iQ ^dF ^ dG_ 1 SF, ^ do 
dt dx dz Re dx dz 


' p' 


pu 


pw 

pu 

F = 

pu^ + p 

G = 

puw 

pw 


puw 


pw + p 

e 


{e + p)u 


(e + p)w 


(2.7) 


and 


F„ = 


'O' 


'O' 



T 

XX 


XZ 



T 

xz 


zz 

M. 


.<?4. 




^zz=jM(>v,-l/2u,) 


=U(w, +uj 

/4=«^xx + W^« + 




g4=«T^+WT^ + 


Re = 


Pr(r-l) 

M 

Pr(r-1)‘ 
UL 


U«, is the free stream velocity and L is the reference length. The pressure is related 
to the other variables by 
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In the above equations the ratio of the specific heats, yis 1.4 and a is the local speed of 
sound. The density is non-dimensionalized by the free stream density , the velocities 
by the free stream speed of sound , and the total energy by . 


E. TURBULENCE MODEL 

The unsteady Navier-Stokes equations can completely model the fluid flow, but 
the computational calculations of turbulent flows for realistic geometries at high 
Reynolds numbers demand very high grid densities and very small time steps. Therefore, 
in order to compute turbulent flows for configurations of practical interest, turbulence 
modeling is used. Turbulence models are implemented with the time-averaged forms of 
the Navier-Stokes equations. 

The most commonly used averaging procedures are: 

• the standard time averaging procedure for incompressible flow and 

• the mass-averaged approach for compressible flows. 

The time averaging procedure omits the high frequency information of the 
turbulence, but the unsteady mean flow information is preserved. For the incompressible 
case the randomly- changing flow variables are replaced with their averages plus their 
fluctuations. Thus the u velocity component is represented as u-u + u' where u is the 
mean velocity and u' is the fluctuation about the mean. The governing equations are time 
averaged and the average of the fluctuation terms is set equal to zero. 

In the compressible flow case the mass-weighted variable of the Favre averaging 

approach is used. In this case u is represented u = u + u" where u is 

_ pu 
u =-^ 

P 

Here time average of the doubly primed fluctuating quantities is not equal to zero. 


( 2 . 8 ) 


After the substitution is carried out for all of the fluctuating flow variables, the entire 
equation is time averaged. Next, all the time averaged terms that are doubly primed 
and multiplied by density are defined to be zero. For example 

pu'' = 0 ( 2 . 9 ) 

The equations of mean motion resulting from the time averaging procedure have 
more unknowns than equations. This constitutes the closure problem of turbulence. 

In order to close these equations a turbulence model must be used. 

The Baldwin-Lomax (B - L) turbulence model [Ref. 10] is a two-layer eddy 
viscosity model which simulates the effect of turbulence in terms of the eddy viscosity 
coefficient |if The term p in the stress terms is replaced by p+pt and the p/Pr in the heat 
flux terms is replaced with pyPr+Pt/Pn- The B-L turbulence model bypasses the need for 
finding the edge of the boundary layer by using vorticity instead of the boundary layer 
thickness. This model is adequate for flows which have mild pressure gradients, but it is 
not very suitable for highly separated flows. The basic equations of the model follow. 

In the inner layer, the eddy viscosity, is assumed to be proportional to the mixing 
length squared and vorticity, and in the outer layer it uses an exponentially decaying 
formula. The irmer eddy viscosity is computed up to the point where it is equal to the 
outer eddy viscosity as shown below. 


(/^I) inner for 

y ^ y crossover I 

) outer for y> 

ycrossover J 

where y is the normal distance from the wall and ycrossover taken at its minimum value 
where it equals y. The inner eddy viscosity is given by: 

W inner 


where 


i = m 


l-exp(-—) 


\coH= 


d(0 du ^ 
V 
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is an experimentally determined damping constant, k is the Von Karman constant. 
The outer eddy viscosity is given by 





III. NUMERICAL FLOW SOLUTIONS 


A. INCOMPRESSIBLE INVISCID FLOW 

The basic governing equation for two-dimensional, incompressible, irrotational 

92^ 9 ^<I) 

inviscid flow is the Laplace equation r- -i- —^ = 0. This comes from the continuity 

dx^ dz 

Eq.(2.2) for incompressible flow (p=const) and the equation for irrotational flow 

du _dw 
dz dx 

A computer code developed by Teng [Ref. 15] was used for the flow calculations. 
This panel code uses the assumption that the airfoil can be approximated by a number of 
panels each having a source strength qj. All the panels are assumed to have the same 
vorticity strength y. These sources and vortices induce a velocity at every point around 
the airfoil. For the flow solution to be found two conditions are applied: 

• The flow tangency condition which states that the normal component of the 
total velocity at the midpoint of each panel due to freestream velocity and the velocities 
induced by the sources and vortices on all the panels is zero. 

• The Kutta condition which states that the pressure on the upper and lower 
surface of the trailing edge must be equal. 

For an oscillating airfoil an additional wake panel is assumed to be attached to the 
trailing edge. The airfoil motion is divided into small steps. At each step the vorticity of 
the wake panel is assumed to be concentrated into a single point vortex which detaches 
from the airfoil with the local velocity such that the Helmholtz theorem can be applied: 
the change in circulation about the airfoil between time steps k-1 and k is equal and 
opposite in direction to the vorticity released into the wake or 
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where 


A is the wake panel length 
Yw is the vorticity strength and 
r is the circulation about the airfoil 

Two more equations are needed in order to determine the length and the 
orientation of the wake panel. Thus two additional conditions must be specified: 

The wake panel is oriented in the direction of the local resultant velocity at its mid-point 
The length of the wake panel is proportional to the magnitude of the local resultant 
velocity at its mid-point and the time-step. Figure 3-1 shows the unsteady wake and its 
elements. 


A. V (Ik-2”Ik-l) 





Figure 3-1 Incompressible Flow wake model 


B. COMPRESSIBLE FLOW 

The governing equations of the compressible flow over airfoils are so complicated 
that no exact solution can be obtained analytically. Their solution is accomplished with 
numerical methods which were applied on a mass scale after the advent of high speed 
computers. 
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The total process of obtaining flow solutions about problems involving fluid 
motion can be represented schematically in Figure 3-2. 


FOR EACH ELEMENT OF FLUID 


Conservation of mass => Continuity Equation 
Newton’s second law => Euler Equations 

Navier-Stokes Equations 
Conservation of energy => Energy Equation 


Solve the equations+B.C’s+LC’s 


1 

r 

Velocity Distribution : 

u(x,z,t), w(x,z,t) 

Pressure : 

p(x,z,t) 

Density : 

p(x,z,t) 

Temperature : 

T(x,z,t) 


r 

Deduce flow behavior 

flow separation 
flow rates 
heat transfer 
forces on bodies 
(skin friction, drag, lift) 


Figure 3-2 Overview of Computational Fluid Dynamics 


The governing partial differential equations, are replaced with systems of 
algebraic equations, so that a computer can be used to obtain the solution. The process of 
converting the continuous governing equations to a system of algebraic equations is 
known as discretization. 
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For finite differencing schemes, a grid of discrete points is distributed throughout 
the computational domain in time and space. The algebraic equations link together values 
of the dependent variables at adjacent grid points. The required number of grid points for 
an accurate solution typically depends on three factors; the dimensionality, the geometric 
complexity and the severity of the gradients of the dependent variables. At each grid 
point each dependent variable and certain auxiliary variables must be stored. All of these 
variables must be stored in main memory for the computation to be efficient. 

As the governing equations for flows around an airfoil are nonlinear, the process 
of the computational solution must be iterative. The iterative process is often equivalent 
to advancing the solution over a small time step. The number of iterations, or time steps, 
might vary from a few hundred to several thousand. 

The discretization process always introduces an error. As long as the discrete 
equations are faithful representations of the governing equations, this error can be 
reduced by refining the grid. If the numerical algorithm that performs the iteration is 
stable, then the computational solution can be made arbitrarily close to the true solution 
of the governing equations by refining the grid. 

1. Grid Generation by Algebraic Mapping 

In order to use an unweighted differencing scheme the flow equations must be 
transformed to a generalized coordinate system using the following transformations for 
the case of unsteady two-dimensional flow: 

4=^(x,z,t) 

C=Cix,z,t) ( 3 . 1 ) 

z=z(x,z,t) 

The above equations can be used to transform the governing equations from the 
physical (x, z,t) to the computational domain (^,^,x) with the Jacobian transformation. 

For more details refer to [Ref.2]. 
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Grid generation has to do with the establishment of the correspondence between 
points (x, z) in the physical domain and points (^,Q in the computational domain. 

Algebraic mapping techniques interpolate the boundary data in one or more 
dimensions in order to generate the interior grid. The generated grid should be well- 
conditioned, i.e., smoothly varying, close to orthogonal and with local grid aspect ratios 
close to unity. For fluid flow problems, the solution is often changing rapidly close to a 
particular surface. It is important to construct a grid that is orthogonal, or near- 
orthogonal, adjacent to such a surface. 

Stretching functions on the boundaries are used in order to define the grid points 
in the interior. Consequently, the two-boundary and multisurface techniques are still able 
to get smoothly viarying near-orthogonal grids with only one-dimensional explicit 
interpolation. 

The distribution of points along the boundary of the domain is handled effectively 
by defining normalized, one-dimensional, stretching functions along boundary segments, 
typically corresponding to each side of the computational rectangle in the (^, 77 ) plane. 
Boundary stretching functions are applicable whether the interior grid is generated by 
solving a partial differential equation or by an algebraic mapping. 

2. Numerical Implementation of Algebraic Mapping 

Computer programs have been developed capable of generating grids between 
two bounding curves based on algebraic mapping techniques. 

For the domain shown in (Figure 3-3), the bounding surfaces consist of a 
symmetric slender body extended downstream, ABC, and a farfield boundary, FED. 
Between these two boundaries half of a C-grid is to be generated. By symmetry the 
complete C-grid can be obtained by reflection about the x-axis. 

As mentioned before, grid generation is split into two parts: 

First, grid point locations on all boundaries are determined and a 1-D stretching 
function is used to control the distribution on the boundaries. 
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Subsequently, the interior grid is generated by the multisurface technique. Two 
intermediate surfaces, Z2, and Z3, are introduced, one each adjacent to the bounding 



surfaces ABC and FED. The. parametric (r) correspondence of surfaces Z 2 , and Z 3 to their 
neighbouring bounding surface is adjusted so that grid lines intersect the bounding 
surfaces orthogonally. The mechanism of choosing x(r), y(r) on surfaces Z 2 , and Z 3 
requires an orthogonal projection, conceptually similar to the near-orthogonal grid 
constmction. 


3. CFD Techniques 

Computational techniques are used to obtain an approximate solution of the 
governing equations and boundary conditions. For example, for two-dimensional 
unsteady incompressible flow, velocity and pressure solutions, u(x, z, t), w(x, z, t) and 
p(x, z,.t), would be computed. The computational solution is obtained in two steps that 
are shown in Figure 3-4. 

In the first step, the partial differential equations and boundary and initial 
conditions are converted into a discrete system of algebraic equations. This step, as 
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mentioned before, is called discretization. The second step comprises the solution of the 
system of algebraic equations. 

The fact that the differentiated terms in the governing partial differential 
equations are replaced by algebraic expressions connecting nodal values on a finite grid 
introduces an error. In order to minimize this error, the most appropriate algebraic 
expressions should be chosen. Equally important as the error in representing the 
differentiated terms in the governing equation is the error in the solution. 



Figure 3-4 Overview of the computational solution procedure 


a. Discretisation process example 

To convert the governing partial differential equation(s) into a system of 
algebraic equations (or ordinary differential equations), a number of choices are 
available. The most common are the finite difference, finite element, finite volume and 
spectral methods. 

The way the discretisation is performed also depends on whether time 
derivatives (in time dependent problems) or equations containing only spatial derivatives 
are being considered. In practice, time derivatives are discretized almost exclusively 
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using the finite difference method. Spatial derivatives are discretized by either the finite 
difference, finite element, finite volume or spectral method. 

The discretization process can be illustrated by considering the scalar 

equation 

du du _ 

—+ «—= 0 (3.2) 

ut ox 

The most direct means of discretization is provided by replacing the 
derivatives by equivalent finite difference expressions. It is based on forward difference 
in time and a central difference in space and gives, after multiplying with At: 

( 3 . 3 ) 

2Aii: 

This equation is referred to as explicit in time. That is if we know u " at all 
grid points in space at time level n, then we have a series of explicit calculations to 
deteimine u""^^ at all grid points. Further, all values of Ui”'*’^ are obtained independently of 
each other; they depend only on Uj", not on the other Ui""^^ terms. 

The process of discretizing Eq.(3.2) to give Eq.(3.3) implies that the 
problem of finding the exact (continuous) solution u(x, t) has been replaced with the 
problem of finding discrete values u”, i.e. the approximate solution at the (i, n)th node. In 
turn, two related errors arise, the truncation error and the solution error. 

The precise value of the approximate solution between the nodal (grid) 
points is not obvious. Intuitively, the solution would be expected to vary smoothly 
between the nodal points. In principle, the solution at some point (Xr, Ur) that does not 
coincide with a node can be obtained by interpolating the surrounding nodal point 
solution. 

To provide the complete numerical solution at time level (n + 1), Eq. (3.3) 
must be applied for all the nodes i=2,...I-l, assuming that Dirichlet boundary conditions 
provide the values ui""^^ and u""^'. 

The discretization process invariably introduces an error unless the 
underlying exact solution has a very elementary analytic form. In general, the error for a 
finite difference representation of a derivative can be obtained by making a Taylor series 
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expansion about the node at which the derivative is being evaluated. The evaluation of 
the leading term in the remainder provides a close approximation to the error if the grid 
size is small. However, the complete evaluation of the terms in the Taylor series relies on 
the exact solution being known. 

b. Convergence 

A solution of the algebraic equations that approximate a given partial 
differential equation is said to be convergent if the approximate solution approaches the 
exact solution of the partial differential equation for each value of the independent 
variable as the grid spacing tends to zero. Thus we require 

u" —^u(Xi,t^) as Ax,Ar—>0 

where u is the exact solution 

The difference between the exact solution of the partial differential 
equation and the exact solution of the system of algebraic equations is the solution error, 
denoted by e; that is 

e," 

The exact solution of the system of algebraic equations is the approximate solution of the 
governing partial differential equation. The exact solution of the system of algebraic 
equations is obtained when no numerical errors of any sorts such as those due to round¬ 
off, are introduced during the computation. The magnitude of the error, at the (i,n)!*' node 
typically depends on the size of the grid spacings. Ax and At, and on the values of the 
higher-order derivatives at that node, omitted from the finite difference approximations to 
the derivatives in the given differential equation. 

4- Navier - Stokes Solution 

If we apply the generalized transformation to the compressible Navier-Stokes 
equations written in vector form Eq.(2.7), the following transformed equation 
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is obtained 


dt'^d^'^dC ReV^ 

where Q is the conservative variables vector: 

P 

pu 
pw 
e 

J is defined in Eq.(3.9) 

F, G are the inviscid flux vectors: 

pu 1 \ pw 

puU + ^^p puW + ^^p 

pwU + ^^p J pwU + ^^p 

{e + p)u-^,p\ \_{e + p)W -C,p_ 

and S is the thin layer approximation of the viscous fluxes in the ^ direction: 

0 

pm^u^ + (p/3)m2C^ 
pm^w^ +(///3)m2^, 
pm^m^ +{pl3)mpn^ 

where 

m,=(u^ + w^)/2+ ^ 

' Pr(r-l) 

m^=Lu + Cz^ 






(3.4) 


The U and W are the contravariant velocity components, tangential to the constant ^ and 
^ surfaces, respectively, and they are given by: 
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u=u^^+w4,+4, 

W = uC,-^w(^+(, 


Non-dimensionalization of the governing Eq.(3.4) is performed using the 
following reference quantities: density p„, length c, time c/ac, and energy po,a^„ 

After non-dimensionalization of the governing equations, Euler solutions are 
obtained setting the viscous terms of Eq.(3.4) zero and applying the flow tangency 
condition at the airfoil surface. For Navier-Stokes solutions the no-slip condition is 
applied at the surface. 

The numerical integration of the Eq (3.4) is performed using an implicit scheme 
[Ref. 11] given by 




+ '>; - Re-' f x fe' - Q', )= 


(3.5) 


■ 2a (^a+1/2 “'^a-1/2 )J 


In Eq.(3.5) V, A and 6 are the. forward, backward and central difference operators 
respectively, whereas the h^=At/A^ , =AtlA^ , A* =dF/dQmd 

B~ =dGldQ are the flux Jacobian matrices. The superscript (.)” refers to the time step 
and the superscript (.)^ refers to Newton subiterations within each time step. 

Time accuracy of the implicit numerical solution is obtained by performing 
Newton iteration to convergence for each time step. Linearization and factorization errors 
are minimized because the left hand side of Eq.(3.5) can be driven to zero at each time 
step. Typically two to three subiterations are sufficient to drop the residuals two orders of 
magnitude during the Newton iteration process. 


19 




,20 



IV. AEROELASTIC ANALYSIS 

A. FLUTTER - GENERAL 

Flutter can be defined as the dynamic instability of an elastic body in an airstream. 
The aircraft aerodynamic surfaces, such as wings, tails, and control surfaces (i.e. bodies 
subjected to large lateral aerod 5 mamic loads of the lift type) are the bodies that are most 
likely to experience flutter. The only air forces necessary to produce it are those due to 
deflections of the elastic structure from the undeformed state. 

The flutter speed Up and frequency (Op are defined, respectively, as the lowest 
airspeed and the corresponding circular frequency at which a given structure flying at 
given atmospheric density and temperature will exhibit sustained, simple harmonic 
oscillations. All small motions are stable at speeds below Up, whereas divergent 
oscillations can occur in a range of speeds (or even at all speeds) above Up. In other 
words, flight at Up represents a borderline condition or neutral stability boundary above 
which flutter will happen. 

. The free vibration of a linear structure in vacuum is a real, or single, eigenvalue 
problem. In contrast, the theoretical flutter analysis leads to a complex, or double, 
eigenvalue problem, where two characteristic numbers determine the critical flutter speed 
and frequency. This is because the simple harmonic motion assumption is made; all 
dependent variables are proportional to e'“‘, and the problem is to find all the 
combinations of U and co for which flutter actually occurs. 

There are three main ways in order to study and predict flutter: 

• theoretical flutter computation 

• wind-tunnel experiments on scaled dynamic models, 

• flight testing of full-scale aircraft. 
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The decision as to which of these is most economical in a given case is very difficult and 
depends on a multitude of factors such as: 

• the purpose and the phase of the study 

• anticipated margin of safety from flutter, 

• the Mach number range, 

• the number of different mass and sfructural configurations to be analyzed. 
Current regulations require that flight flutter testing must be preceded by flutter analysis 
and wind tunnel testing. 


B. THE PHENOMENON OF FLUTTER 

The phenomenon of flutter can be described with a simple experiment: consider a 
cantilever wing, mounted in a wind tunnel at a small angle of attack and with the root 
rigidly built in. Then when there is no flow in the tunnel we deliberately deflect and then 
release the wing. We observe that the oscillation decays gradually. Then we increase the 
speed of flow in the wind tunnel and repeat the experiment; the oscillation of the 
disturbed airfoil dies out. With further increase of the flow speed, however, a point is 
reached at which the oscillation maintains itself with steady amplitude. This is the critical 
flutter speed. If the wind tunnel speed is increased, flutter will occur with increasing 
amplitude. It is obvious that during the aircraft flight at speeds above the critical, an 
initial disturbance (coming from a gust for example) will cause an oscillation of great 
violence. In such circumstances the airfoil suffers from oscillatory instability and is said 
to flutter. 

The above experiment on wing flutter shows that the oscillation is self-excited; 
i.e., no further external disturbance is required. The motion can maintain itself or grow 
for a range of wind speeds, which depends on the design of the wing and the conditions 
of the test. 
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An aerodynamic surface, as an elastic body, has infinitely many degrees of 
freedom (d.o.f). But often it is sufficient to consider only three variables, namely the 
flexure (plunge), the torsion (pitch), and the control-surface rotation. A flutter mode 
consisting of all three elements is called a ternary flutter. In special cases, however, two 
of the variables predominate, and the corresponding flutter modes are called binary flutter 
modes. In fact, many airplanes can be replaced by a system of simple beams, so that the 
elastic deformation can be described by the deflection and torsion of the beams, in 
addition to the rotation of control surfaces about their hinge lines. 

For the oscillatory motion of a wing only the first two components of flutter are 
often considered: the flexural (plunge) and the torsional (pitch). In general, the coupling 
of those two degrees of freedom is an essential feature for wing flutter. The oscillation 
that occurs at the critical speed is harmonic. Experiments on cantilever wings show that 
the flexural movements at all points across the span are approximately in phase with one 
another, and likewise the torsional movements are all approximately in phase, but the 
flexure is considerably out of phase from the torsional movement. As will be seen later, it 
is this phase difference that is responsible for the occurrence of flutter. 

A rigid airfoil in a low-speed attached flow which is constrained in torsion and 
can only plunge (flexural d.o.f) does not flutter. In contrast, a rigid airfoil with only the 
torsional d.o.f can flutter for some special mass distributions and elastic-axis locations. In 
this text we will restrict the term “flutter” to the oscillatory instability in a flow which 
does not experience flow separation or strong shocks. 

For a control surface, such as a flap or an aileron, its flexural and torsional elastic 
deformation are not so important as its freedom to turn about the hinge line. 

Consequently, the deflection of a control surface can simply be described by the angle of 
rotation about its hinge line. 

The degrees of freedom of flutter mentioned before, together with the freedom of 
the airplane to move as a rigid body, offer a large number of possible combinations of 
binary, ternary, and higher modes of flutter. Since it is not clear which of these modes 
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correspond to the actual critical speeds, it is necessary either to resort to experiments or, 
in a theoretical approach, to analyze all cases. This is why a successful flutter analysis 
depends so much on the analyst’s experience. He must be able to choose, among all 
possible modes, those that are likely to be critical for a given structure. 

As one can easily understand, flutter analysis is an extensive subject; for this 
reason it is preferable to divide its presentation in sections. Some basic definitions are 
given in Section C. The role of the elastic stiffness in flutter prevention are explained on 
an empirical basis in Section D. The origin of flutter from the aerodynamic point of view 
is then considered in Section E. It will be shown that flutter occurs because the speed of 
flow affects the amplitude ratios and phase shifts between motion in various degrees of 
freedom in such a way that energy can be absorbed by the airfoil from the airstream 
passing by. The physical explanation and the equations of motion for the one d.o.f flutter 
will be discussed in Section F, whereas the two d.o.f flutter will be explained in Section 
G. A practical way to find the flutter condition with the Panel Code (UPOT) will be 
shown in Section H. 


C. NONDIMENSIONAL PARAMETERS 

In flutter analysis any non-dimensional quantity relating to the motion can be 
expressed as a function of two parameters: • 

p/(j and K/dfU^ 

where 

1 Typical linear dimension 

U Air speed 

p Air density 

a Typical density of structural material 

K Typical torsional stiffness constant (ft-lb per rad) 
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If we consider the energy dissipation of the structure, the viscosity and the 
compressibility of the fluid in flutter analysis then three more non-dimensional variables 
must be added: 

g Material damping coefficient 

Re Reynolds number 

M Mach number 

For two systems to be dynamically similar they should have equal variables g. Re, 
and M in addition to the above parameters. In general, g is important in control-surface 
flutter. Re is important in stall flutter, and M is important in high-speed flight; otherwise 
their effects are small. 

The non-dimensional form of the frequency of oscillation co (radians per second, 
with dimension T^) is: 



which is called the reduced frequency. 

Two similar systems having the same values of yo/crand K/al^U^ flutter at the 
same reduced fi-equency. Since all derived concepts relating to the motion can be 
expressed in functional relations as above, it is clear that the equality of the values of the 
parameters /?/(Tand K/al^U^ is sufficient to guarantee dynamic similarity of the two 
systems. 

The reduced frequency characterizes the variation of the flow with time. Its 
inverse, U! al is called the reduced speed. The physical meaning of the reduced 
frequency is given as follows: Consider that a disturbance occurs at a point on a body and 
oscillates together with the body. The fluid element influenced by the disturbance moves 
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downstream with a mean velocity U. Let the frequency of oscillation of the body and the 
disturbance be (O. Then the spacing, or "wave length" of the disturbance, is'k=T-U or 
X=2nU/(0. So the reduced frequency for the airfoil of chord c is: 


, CQC . C 

k- — = 2n:— 
U A 


(4.2) 


Therefore, the reduced frequency represents a ratio of the characteristic length of the body 
(chord) to the wave length of the disturbance. In other words, the reduced frequency 
characterizes the way a disturbance is felt at other points of the body. Since every point 
of an oscillating body disturbs the flow, one may say that the reduced frequency 
characterizes the mutual influence between the motion at various points of the body. 

The reduced natural pitching frequency is defined as 

(o„c 


k„ =■ 


U 


(4.3) 


where 


K. 


0>a=. 


(4.4) 


is the uncoupled natural torsional frequency of the system, K„is the spring constant for 
pitching and I„ is the moment of inertia about the pivot point per unit span. From Eq.(4.3) 
and (4.4) we have: 



(4.5) 


The reduced flutter velocity is defined: 



U 

(0„c 



(4.6) 
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D. STIFFNESS CRITERIA 


As mentioned before, wing-flutter occurs when the reduced frequency is lower 
than the following critical value: 


k 


cr 


U 


where U is the mean speed of flow, co„ is the undamped natural pitching frequency of the 
wing, and c is the chord length of the vibrating portion of the wing. 

For safety against flutter, the reduced frequency should be higher than kcr, or the 
design speed of the airplane should be lower than 




(4.7) 


The frequency co^, which may be determined by ground vibration experiments or 
computed by a theoretical analysis, increases with increasing stiffness of the stracture. 
Therefore, the critical speed can be raised by increasing the wing stiffness. It is obvious 
that when the lower limit of Ucr is defined (e.g., the maximum speed of flight) we can 
determine the minimum value of co, and hence the minimum value of the torsional 
rigidity. 


E. VERTICAL TRANSLATORY OSCILLATION 

The fundamental role played by the aerodynamic forces in inducing flutter is that 
they are the means through which kinetic energy is transferred from the airstream to the 
airfoil as elastic and kinetic energy. Hence, the possibility of flutter can be discussed by 
considering the energy relation. From the analysis we can find the amount of energy 
exchanged between the airstream and the airfoil; if the oscillating body gains energy from 
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the airstream in completing a cycle then the oscillation is aerodynamically unstable; if it 
gives energy to the airstream the oscillation is going to be aerodynamically stable. 

In order to determine the energy transfer consider an airfoil performing a vertical 
translatory oscillation with a constant amplitude ho. Then the vertical displacement can 
be described by 

h = h,e“^‘ (4.8) 

where h is defined to be as positive downward. The speed of downward motion is 
therefore 

h = icohQe^"' ' (4.9) 

where the prime indicates a differentiation with respect to time. If A is a constant, the 
downward motion will induce a lift force Lq on the airfoil: 


Lo=\pU^S 


dCj^ h 
~daU 


(4.10) 


where _is the dynamic pressure and S is the wing area. This value Lg is the quasi¬ 

steady lift and is defined as positive upward in the usual sense. As the airfoil performs a 
translatory oscillation, the true instantaneous lift acting bn the airfoil differs from Ig both 
in magnitude and in phase and its value can be stated in the form 

L = V^'^ (4.11) 

In this expression r represents the ratio of the absolute value of the instantaneous lift to 
that of the quasi-steady lift, and \|; the phase angle by which the actual lift leads the quasi¬ 
steady value. The quantities r and X)/ depend on the reduced frequency k, the Mach 
number M, and the Reynolds number Re. For a nonviscous incompressible fluid, r and \|r 
are functions of k alone. The ratio L/Lg = re‘'^ can be plotted as vector with length r and 
angle \|/. Using incompressible oscillatory thin airfoil theory von Karman and Sears 
derived the vector diagram of Figure 4-1. 
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Figure 4-1 Vector diagram of lift in vertical translation oscillation. (From Y.C.Fung Ref. 1) 


When the airfoil moves through a distance dh, the differential work done by the 

lift is 


dW = -Ldh = -Lhdt (4.12) 

and therefore the total work done by the air on the airfoil during one cycle of oscillation, 
as given by Y.C Fung, [Ref. 1] is: 


ItcIco j 27t(a) 

W = - f RG[L]-K&[h]-dt = -—pUS—^{(johQfr [ sin{oX + \i/)-sm{(jOt)dt = 
n 2 da i 


K dC, 9 

= — pUS — —cch^ rcosy/ 
2 da 


(4.13) 


Hence, the gain of energy Wby the airfoil from the airstream is proportional to 
-cost]/, where \j/ is the phase angle by which the actual lift leads the quasi-steady value. If 
-7!/2 <\jt < nfl, then Wis negative; i.e., the oscillating airfoil will give energy to the 
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airstream and the oscillation is therefore stable. If we refer back to Figure 4-1, it is seen 
that the condition -n/2 <\)/ < 7i/2 is satisfied. Hence, in a nonviscous incompressible 
fluid, the vertical translation oscillation is aerodynamically stable. 

This example shows the importance of the phase angle between the aerodynamic 
force and the oscillatory motion. Based on this phase angle we can conclude that purely 
translational flutter in one degree of freedom is impossible. 

However as shown in the next Section, this conclusion does not hold for single- 
degree-of-freedom torsional motion. 

F. ONE DEGREE OF FREEDOM FLUTTER (TORSIONAL) 

In this section one degree of freedom torsional flutter will be examined. In 
Subsection 1 the negative damping at low reduced frequency k will be explained 
qualitatively. In Subsection 2 the time and frequency domain analysis will be presented. 

1. Physical Explanation 

• Consider an airfoil in an incompressible airflow of speed U, subjected to a slow 
(low K), nearly sinusoidal torsional motion. Assume that the rotational axis is at the 
leading edge (L.E.), but similar results are obtained if it is anywhere in front of the quarter 
chord point. 

Examine now step by step how flutter may occur for this one degree of freedom 
oscillating airfoil with the help of Figure '' given by Bisplinghoffet al [Ref.4]. At stage 
(a) when the airfoil is at the highest angle oi attack all the vortices shed during the 
preceding last half cycle have counter-clockwise vorticity. Between (a) and (c) the airfoil 
pitches down resulting in vortices of clockwise vorticity. As k is low, the counter¬ 
clockwise vortices are far away from the trailing edge (T.E.) and do not contribute 
signific^tly in inducing the upwash. Therefore, as the airfoil pitches through the zero 
AOA position, stage (b), a lift L 2 is induced which is in the same direction as the motion 
and therefore net positive work is done on the airfoil. 
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When the airfoil increases the angle of attack from steps (c) to (a) the counter¬ 
clockwise vortices dominate inducing a wake upwash and therefore a lift L 2 and moment 
which is in the same sense as the angular velocity and again does positive work on the 
wing. 

This transfer of energy from the airstream to the airfoil is the origin of the 
negative damping. At higher values of k the wave length of the vertical wake becomes 
shorter so that vortices from previous half-cycles can contribute. The net effect is an 
eventual change of the aerodynamic damping from negative to positive. 
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The phenomenon of upwash and lift induced by the wake vortices can be 
visualized by the incompressible panel code UPOT. The airfoil wake vortices are shown 
in Figure 4-3 as squares. The solid ones are the clockwise vortices whereas the dotted 
ones are the counterclockwise vortices. The airfoil is shown in the same four positions as 
in Figure 4-2. It is clearly seen that the rotation of the vortices sketched by Bisplinghoff et 
al is reproduced by the panel code. 



Figure 4-3 Incompressible flow wake visualization (from UPOT) 
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The essence of the foregoing discussion is the inclusion of unsteady effedts. If 
only quasi-steady aerodynamics had been used (without vortex shedding), the possibility 
of torsional flutter would have been overlooked. For a more detailed analysis refer to 
Bisplinghoff et al [Ref. 4 pages 262-263]. 

2. Flutter Analysis 

Flutter analysis can be done with two methods. The first is the time-domain 
analysis in which the airfoil motion is found as a function of time. It is ba^ed on 
discretisation of the ordinary differential equation of motion. The second method is the 
classical theoretical approach based on frequency-domain analysis with which the critical 
flutter speed and frequency are calculated. For both methods we consider a simple system 
of a two-dimensional, rigid airfoil of unit span, which is hinged at a specific point (pivot) 
but elastically restrained from rotating about this axis by a torsion spring with constant 
K„. The airfoil, which has a moment of inertia about the pivot point I„, is placed in an 
airstream so that the unstrained position of the spring corresponds to zero angle of attack 
a (Figure 4-4). Then we let the airfoil perform a torsional oscillation. 
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a — s 


( 4 . 18 ) 



(4.19) 


This is the system of equations used for the incompressible flutter analysis 
where Cm is computed at each time step using the panel code UPOT. 

For the compressible flutter analysis the time can be non-dimensionalized 
using the free stream speed of sound, , instead of the freestream velocity U„. Thus the 
dimensionless time can be defined as 


T = — (4.20) 

c 

Using this non-dimensional Cm, i„and x, the equation of motion (4.19) becomes: 

s' = -kla.+—Mlc„ 
ijr 


The reduced natural frequency k„for this case is based on the freestream 
speed of sound, whereas in the incompressible case it is based on the freestream velocity. 
In order to have compatibility between the values of k„ found for incompressible and 
compressible flow we have to convert the latter into the former values. The procedure 
will be shown in Section H. 

Integration of Eq.(4.18) and (4.19) can be done with a first-order Euler 
scheme, a second order modified Euler scheme or a fourth order Runge-Kutta scheme. 

The second-order modified Euler scheme is derived by applying the 
trapezoidal rule to integrate an equation y = f{x,t) as follows: 

=}'„ +^/2l/(y„^j,t„^,) + /(y„,rj] 
and using an initial value for y(0). 

The fourth-order Runge-Kutta scheme is derived by applying the following 

equation: 

=}'„+(1/6)[A:, +2^2+2^3-1-^:4] 

with 
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K =¥iyn^^n) ^2 =¥(yn + / 2,t „ + k / 2) 

k^=¥iyn + k2^2,t„+h/2) k4=hf(y„+kj,t„+h) 

Panel codes, like UPOT, typically use large time steps and therefore use of 
the fourth-order Runge-Kutta scheme is advisable. In contrast, the stability limits of 
Euler/N-S codes typically require small time steps and therefore the first-order Euler 
scheme gives satisfactory accuracy. 


b. Frequency domain analysis 


In order to perform the classical flutter analysis in the frequency domain 
we assume that the airfoil executes a torsional oscillation with a constant amplitude Oo. 
The angle of attack a at any time can be described by 

a = (4.21) 

and after substituting (4.16) into (4.14), and dividing through by 7 rpb*CO^aQe'‘“we. get 

v2^ 


Ttpb 




+ C„ =0 


(4.22) 


where is the undamped natural frequency of torsional vibration and Cm is the 
dimensionless aerodynamic coefficient which is a complex number. 

Hence (4.18) can be split into real and imaginary parts 

L 


Re{c„> 


Tupb 




0 ) 




(4.23) 

(4.24) 


Flutter occurs only at those values of the reduced frequency that make the 
out-of-phase component of the aerodynamic moment Eq.(4.24) vanish, provided that the 
corresponding in-phase part is of such magnitude that Eq.(4.23) yields a non-imaginary 
flutter frequency (0. The latter condition is met when the rotational axis is ahead of the 
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one quarter chord line. Hence by this method we can find the flutter airspeed and 
frequency. 

The variation of the real and imaginary parts of the aerodynamic moment, 
based on Theodorsen’s thin airfoil theory, for k defined as coc/2U_, is shown in Figure 4- 
5. For the definition ok k according to Eq.(4.2), the imaginary moment becomes zero at a 
reduced frequency of k=0.076. This is consistent with the previous discussion which 
showed that torsional flutter occurs only at low reduced frequencies. For further details 



Figure 4-5 Variation with reduced frequency k of the real and imaginary parts of the dimensionless 
. aerodynamic moment nia due to pitching of an airfoil about its leading edge in 

incompressible flow (for k defined as ox:/2U«). 

G. TWO DEGREE OF FREEDOM FLUTTER 
1. Physical Explanation 

For the two-degrees-of-freedom flutter we will consider two cases: a) the torsional 
deflection leads the bending deflection by 90 degrees as shown in Figure 4-6 whereas b) 
the torsional deflection and bending are in phase with each other as shown in Figure 4-8. 
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In both figures the solid arrows indicate the airfoil motion whereas the dotted arrows 
indicate the airfoil lift. 


For simplicity we will assume purely steady state values for the lift corresponding 
to the instantaneous position of the airfoil. 



Figure 4-6 Two degrees of freedom airfoQ motion - phase difference 90° 









From Figure 4-7 we can see that for the two first quarters of the cycle, both lift L 
and displacement dh are positive whereas for the two last quarters of the cycle, both lift L 
and displacement dh are negative. As a result over the whole cycle work is done on the 
airfoil. 

When the torsional deflection is in phase with the bending deflection then from 



Figure 4-8 Two degrees of freedom airfoil motion - phase difference 0° 




Figure 4-9, we can see that for the first quarter of the cycle, the lift L and the 
displacement dh are positive so positive work is done on the airfoil. For the second 
quarter of the cycle, the lift L is positive whereas the displacement dh is negative so 
negative work is done on the airfoil. Similarly positive and negative work are done on the 
airfoil during the third and fourth quarter of the cycle. As a result zero work is done on 
the airfoil. 

2. Equations of motion 

The equations of motion for the two-degree-of-freedom system with no 
mechanical damping shown, for example, by Fung [Ref. 1] pages 210-211, can be written 
as: 

+ S^a + hK^ = -L (4.25) 

and Sah + Iaa + o^^=M^ (4.26) 

where is the static moment about the elastic axis. 

/„is the airfoil moment of inertia 
L is the airfoil aerodynamic lift 
and is the aerodynamic moment about the elastic axis 

a. Time domain analysis 

For the incompressible flow case the equations of motion (4.25) and (4.26) 
can be non-dimensionalized using the same non-dimensional coefficients that were used 
for the one degree of fireedom case Eq.(4.15)-(4.17). 

Then the system of equations (4.25) and (4.26) becomes: 
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and 


mh" + S^a" + mklh = -2c, Ik 
sy + i^a'' + i„kl=2c^/K 


where the primes refer to differentiation with respect to nondimensional time x. 
. In a matrix form the above system of equations becomes: 


with 


[M][Y}"+H3r}={F} 








(4.27) 


For the compressible flow case (N-S Code) the equations of motion (4.25) 
and (4.26) can be non-dimensionalized using the free-stream speed of sound, , instead 
of the free-stream velocity t/„ such that the dimensionless time can be defined from 
Eq.(4.20). 

Then the system of equations have the same form as Eq.(4.27) with the 
difference that the vector {p} now is defined as 

Eq.(4.34) can be rewritten as: 

{X } = [M ]-' }- [M ]-■ \k Kx } (4.28) 
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Defining {X } = {f}, we can solve Eq.(4.28) as a system of two first-order 
differential equations 

As we mentioned for the one degree of freedom case, the integration of this system of 
equations is done with the fourth-order Runge-Kutta scheme in the panel code, and with 
a first-order Euler scheme in the Euler/N-S codes. 

b. Frequency domain analysis 

Assuming constant amplitude time harmonic oscillation including plunge 

and torsion 

h = 

a = a 

the equations of motion become: 

e“^‘\rCO^Mh,-(0^a,S^+h,K,]=L (4.29) 

e““‘\rO)^I,^oc,-co\S,+aoKj=M, (4.30) 

The aerodynamic lift L and moment are assumed to be linear functions of ho and oco so 
Eqs. (4.30) and (4.31) constitute a system of homogeneous equations. Therefore the 
determinant of these equations must be set to zero in order to obtain a solution. As 
discussed in [Ref 1 and 6] the solution of the flutter determinant normally requires an 
iterative procedure. For further details refer to these text books. 


42 



H. NS - UPOT REDUCED FREQUENCY COMPATIBILITY 


As mentioned before, the reduced frequencies calculated with the Euler and N.S 
codes need to be redefined in order to be compatible with the conventional definition of 
k. This is done as follows: 

The non-dimensional time used in UPOT is defined as 

tu^ 

'^vpor ~ (4.31) 

c 

whereas in the N.S code it is defined as 



(4.32) 


c 

In order to have comparable results from the compressible and incompressible flow 
solutions we have to adjust the reduced frequency values of the NS code to those of the 
UPOT code. From the definition of the reduced frequency: 



(4.33) 


we get 


k — 


In c 

TIT 


where T is the period of airfoil oscillation. 

For the results computed with the UPOT code and because of (4.1) we have 

^ 2n c In 

^vpoT^ Typ^yj. 


Similarly for the NS results from (4.4) and because of (4.2) we get 


(4.34) 


(4.35) 


2n c _ 2n c _2n 1 
T^sC Af„ 


(4.36) 
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The last equation shows that we have to divide the reduced frequencies that come from 
the NS code {In! ) with the Mach number () in each case in order to find the 
reduced frequency which is compatible with the UPOT results. 

I. FLUTTER PREDICTION WITH UPOT CODE 

One practical way to find the flutter condition for a given airfoil with the UPOT 
code is to see the resulting diagram of the moment coefficient with respect to AOA. The 
work is proportional to the integrated area in the loop. Positive area means that work goes 
into the flow. Negative area means that work goes into the airfoil. Zero area means that 
the airfoil flutters. Running the code for low values of reduced frequency it is seen that 
the aerodynamic moment coefficient forms a clockwise elliptic path as the AOA changes 
during the airfoil oscillation (Figure 4-10), so net positive work is done on the airfoil. 



Figure 4-10 Aerodynamic moment coefficient graph for low k values 
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For high values of reduced frequency the path has a counterclockwise direction as 
shown in Figure 4-11 and the work is negative. 



Figure 4-11 Aerodynamic moment coefficient graph for high k values 

At last, when the critical value of reduced frequency is used the path becomes a 
straight line as shown in Figure 4-12, so no work is done on the airfoil. 
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Figure 4-12 Aerodynamic moment coefficient graph for the critical flutter k value 
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V. RESULTS 


A. GENERAL 

All the results presented in this thesis were obtained with the time domain method 
described in Chapter IV, and the flow field solutions at each time step were computed 
based on the Navier-Stokes, Euler or unsteady potential flow (UPOT) solution, described 
in Chapter m. 

The reduced frequency values used in the Navier-Stokes and Euler codes were 
redefined in order to be compatible with the UPOT definition. The grid generation for the 
airfoils and the input data for the Navier - Stokes code are presented in Sections B and C. 
The procedure for extracting the results is described briefly in Section D. Then the 
calculation for the non-viscous (Euler) case is presented in Section E. The calculation for 
a wide range of pivot points is shown in Section F. The effect of Mach number and airfoil 
thickness on torsional flutter is presented in Sections G and H. The influence of the 
viscous effects and a comparison between viscous and non-viscous flow are shown in 
Sections I and J. In Section K a comparison between UPOT and NS results for low 
subsonic Mach numbers (M=0.1-0.3) is presented. 

B. GRID GENERATION 

Airfoil grids for NACA 0006, 0009,0012 and 0015 were computed with 
ALGEM. In Figure 5-1 the NACA 0015 grid for the Euler calculations is presented 
whereas in Figure 5-2 the grid details close to the airfoil are shown. This NACA 0015 
coarse grid size is 201X41 points. 

For viscous flow different grids were used because of the.necessity to resolve the 
flow close to the airfoil. The grid size is 201X61 points. 
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Figure 5-1 C-grid for NACA 0015 (non-viscous flow - every other grid line shown) 


48 


































Figure 5-2 NACA 0015 C-grid details (non-viscous flow) 


C. INPUT DATA DESCRIPTION 

The input data file includes many variables, some of which have to be changed 
during the calculation procedure. A typical input file for the initial calculation of the 
steady solution for the non-viscous flow is presented next. The use of the input file 
variables is also explained in this section. 





INPUT FILE: 

# IREAD ITER NPRESTT NLOAD ODVAR 

0 10000 1 1 1.00 

# ALPHA OSCIL RAMP REDFRE ALFAMND ALFAMXD 

0.5 false false 0.100 -.500 0.500 

# PLUNGE PLMX PLMY PLPHSXD PLFREQ 

false 0. 0.10 0. 0.1 

# MACH RE Vise TURBL 

0.700 0. false false 

# TIMEAC COUR NEWHT 

false 30.00 1 

# free mass ialpha ka kh xp xa 

false 0.0 100.0 0.200 0.0 0.0 0.0 

# aneut hneut hO 

0.0 0.0 0.0 

INPUT VARIABLES EXPLANATION: 

IREAD 0 No initial solution, free stream conditions initialize the 
flowfield in the startup steady-state computations 
1 Initial solution is read from a binary file saved from 

the previous run (default, at the end of each run, solution 
file is saved as binary) 

2 Initial solution is read as formatted (plotSd form) 

-1 Initial solution is read as binary, unsteady motion starts 

-2 Initial solution is read as formatted (plotSd form), , 

unsteady motion starts 

ITER # of timesteps 

NPRINT Residuals are printed out at every nprint timesteps 
NLOAD Aerodynamic loads are printed out in ns.out and written into 
lo.d file at every nload timesteps 

ODVAR Solution variables, q array, are written into "qp.d” file 
at every delta odvar change in unsteady motions 
(degrees in oscillatory motion, amplitude change in plunge) 
ALPHA Steady state AO A (do not set it to zero, instead set it to 0.0001) 
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OSCIL 


false N/A 

trae sinusoidal oscillations in pitch 


RAMP false N/A 

true straight ramp motion in pitch 

REDFRE Reduced frequency of the unsteady pitching motion, 

based on the half chord, chord length is assumed to be 1. 

ALFAMND Min AOA of the pitching motion 


ALFAMXD 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 

REYNOLD 

MACH 

Vise 

TURBL 

TIMEACC 


COUR 


Max AOA of the pitching motion 
false N/A 

true plunging motion 
Plunge amplitude in x 
Plunge amplitude in y 

The phase angle between x and y amplitudes, in degrees 

The reduced frequency of the plunging motion, 
based on the half chord, chord length is assumed to be 1. 

Reynolds number of the freestream flowfield 

Mach number of the freestream flowfield 

false Euler solution 

true Viscous Navier-Stokes solution 

false Laminar flow is assumed 

true Baldwin-Lomax turbulence model is applied 

false Variable local time stepping in the computational grid, used 
in the computations of steady state and/or attached fiowfields. 
true Constant time stepping everywhere in the computational grid, 
used in the computations of unsteady and separated fiowfields 

Courant number of the timestepping (50-1500), determines the 
time step of the computations based on the minimum grid size and 
the freestream conditions. Its value depends on the computational grid. 
For diverging computations, its value needs to be reduced. 
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If the residuals in the output file increases in time, it is the sign that 
Courant number is to be reduced. 

NEWTIT Number of Newton subiterations in each timestep, applied in 
unsteady flows (2-3), for steady flowfields it is set to 1. 

D. CALCULATION PROCEDURE 

For every flutter calculation initially the steady solution had to be found. For the 
non-viscous cases the code was run for 6000 iterations (1000 with COUR Number of 1, 
1000 with COUR Number of 2,2000 with COUR Number of 4 and 1000 with COUR 
Number of 8). For the viscous case and due to the increased number of nodes of the grid 
the calculation needed over 12000 iterations in order to converge to a steady solution. 

The steady solution was assumed to be satisfactory when the step variation of Cl 
coefficient was found to be less than 10'^ or when the five first digits of Cl remained 
constant. 

After the steady solution was found the oscillating airfoil calculation started. 

From Chapter IH Section G we saw that for pivot points behind the quarter chord point 
(0.25c) a statically stable solution is always obtained. Therefore the calculations were 
performed for pivot points ahead of the quarter chord point. The pivot points that were 
chosen for the flutter calculations were: 0.15-c, 0, -0.25-c, -0.50-c, -0.75-c and -0.85-c. For 
e ery case an initial value of the uncoupled reduced natural pitching frequency ka (Eq. 
4.3) was assumed (which corresponds in Eq.(4.4) to an airfoil spring constant Ko) and the 
time variation of the airfoil angle of attack (AOA) was computed. If the AOA was found 
to decrease, the assumed value of k* was adjusted until the motion started to diverge. For 
every computation the slope of the AOA amplitude was determined and an interpolation 
procedure was implemented in order to find the value of k* which gives the critical flutter 
condition. After the critical condition was found the reduced firequency of oscillation was 
calculated. 
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E. NON - VISCOUS (EULER) CALCULATION 


The whole procedure described above will be shown for the non-viscous flow 
over an oscillating airfoil NACA 0015 at M=0.7. The elastic axis of oscillation is the L.E 
of the airfoil (pivot point Xp=0). An initial AOA disturbance of 0.5° will be imposed on 
the airfoil in order to start the motion. 

The input file for the steady solution is as follows (4* run for COUR Number :8): 

# IREAD ITER NPRINT NLOAD ODVAR 

1 2000 1 1 1.00 

# ALPHA OSCIL RAMP REDFRE ALFAMND ALFAMXD 
0.5 false false 0.100 -.500 0.500 

# PLUNGE PLMX PLMY PLPHSXD PLFREQ 

false 0. 0.10 0. 0.1 0.0 0.0 

# MACH RE Vise TURBL 

0.700 0. false false 

# TIMEAC COUR NEWTIT 

false 8.00 1 

# free mass ialpha ka kh xp xa 

false 0.0 100.0 0.00000 0.0 0.0 0.0 

# aneut hneut hO 

0.0 0.0 0.0 

The Steady solution is presented in Figure 5-3. It is seen that Cl reaches a steady value of 
0.1017 after T =5 sec. 


The final result for the steady solution is presented in Table 5-1. 


T 

AOA 

Cl 

Cd 

Cm 


0.5 

0.1017382 

0.0040226 

-0.0257816 


Table 5-1 Steady state solution for NACA 0015 M=0.7 (Euler case) 
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Figure 5-3 Steady state solution for NACA 0015 M=0,7 (Euler case) 

After the steady solution is calculated, the critical value of the reduced natural 
pitching frequency ka for flutter is found, by first assuming an initial value of 0.36, as 
shown in the input file 


# 

IREAD 

ITER 

NPRINT 

NLOAD 

ODVAR 



-1 

10000 

1 

1 

1.00 


# 

ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 


IT) 

O 

false 

false 

0.100 

-.500 

0.500 

# 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 



false 

0. 

o 

H 

o 

0. 

0.1 

o 

o 

# 

DIACH ■ 

RE 

• Vise 

TURBL 




0.700 

0. 

false 

false 



# 

TIMEAC 

COUR 

NEWTIT 





true 

15.00 

3 
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# 

free 

mass 

ialpha 

ka kh 

xp 

xa 


true 

o 

o 

100.0 

0.360000 0.0 

o 

o 

o 

o 

# 

aneut 

hneut 

hO 





o 

o 

o 

o 

o 

o 





As one can see, seven variables changed values comparing with the steady case 
input file: the IREAD was turned to -1, the number of iterations to 10000, the TIMEAC to 
true, the COUR number to 15.0, the NEWTIT to 3 and the free to true. Also ialpha (non- 
dimensional moment of inertia) was chosen to be 100 and the ka to be 0.36. 

The plot of the resulting airfoil motion is presented in Figure 5-4 . 



From the plot we can see that the airfoil motion decays with time so the assumed 
value of ka was too high. In order to determine the critical ka value the maximum AOA 
values of the graph are input into a regression equation to find the slope of the maxitnnm 
AOA values, namely -2.860E-04 (Figure 5-5): 
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Figure 5-5 Time rate of change of AOA amplitude of airfoil motion Ka=0.36 (Euler case) 


Repeating this process for ka of 0.25 we get the AOA history presented in Figure 








From this plot we observe that the airfoil motion diverges. Therefore the assumed 
value of ka is too low. Then applying again a regression equation for the maximum AOA 
values of the graph yields a slope of 4.566E-04 ( Figure 5-7): 



Figure 5-7 Time rate of change of AOA amplitude of airfoil motion Ka=0.25 (Euler case) 

Now we can interpolate the values of the slope found for the two assumed values 
of the reduced natural pitching frequency ka- In order to have more accurate results we 
can follow the same procedure for another value of ka (here assumed 0.30). Finally we 
get the results shown in Table 5-2. 


Ka 

Slope 1 

0.36 

-2.860E-04 

0.30 

1.669E-04 

0.25 

4.5664E-04 | 


Table 5-2 Regression equation slope values for various ka values 


Making a curvefit in EXCEL for the above results (Figure 5-8) we can find that 
the value of the ka for which the slope of the extreme points is zero, is ka= 0 . 324 . 
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y = -0.0068X + 0.0022 



Figure 5-8 Curvefit of vs time rate of change of AOA amplitude (Euler case) 

To verify whether this ka value is indeed the critical reduced natural pitching 
frequency, the code was run for k® =0.324 which yielded the airfoil motion presented in 
Figure 5-9: 



Figure 5-9 Unsteady case solution - Ka=0.324 Critical flutter case (Euler) 
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Omitting the first peaks, which are caused by the starting transient, and applying 
again a regression equation for the maximum points, we find that the slope is 
6.63859M0'^, which is very close to zero (Figure 5-10). Therefore the reduced natural 
pitching frequency ka of 0.324 is accepted as the critical value. If this value didn’t give 
the flutter condition (slope not close to zero), the new slope value would be interpolated 
again until finding the ka for zero slope. 
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Figure 5-10 Time rate of change of AOA amplitude of airfoU motion Kct=0.324 Critical flutter 
. (Euler case) 

The reduced frequency of oscillation is found using (4.36) 

■_ 111 1 

In order to find the period of oscillation the initial transient exhibiting variable 
peaks are ignored, and the average period between the minimum and maximun peaks is 
computed as shown in Table 5-3. 
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Min peaks 

T(NS) 

Max peaks 

T(NS) 

92.6997 


105.954 


119.209 

26.5092 

132.465 

26.51039 

145.721 

26.51206 

158.977 















Average T (NS) = 

26.511 



Table 5-3 Average period T calculation for NACA 0015 flutter - Ka=0.324 


The average period of oscillation (non-dimensionalized) is 26.511. Then applying 
Eq.(4.36) for Mach number of 0.7 we find that the reduced frequency of oscillation is 
0.339. 


The pressure distribution around the airfoil for a whole cycle of oscillation at 30° 
increments is shown in Figure 5-11. 
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AOA=150' 


AOA=180“ 

Figure 5-11 Pressure Contours around NACA 0015 for a whole cycle of oscillation M=0.7 Xp=0.0 
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AOA=330' 


AOA=360' 


Fig.5-11 Pressure Contours around NACA 0015 for a whole cycle of oscillation M=0.7 Xp=0.0 (cont) 
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F. EFFECT OF PIVOT POINT ON TORSIONAL FLUTTER 


The previously described procedure was applied to obtain the effect of pivot 
location, Mach number and airfoil thickness on single degree of freedom torsional flutter. 
Figure 5-12 shows the variation of the flutter reduced frequency with pivot point for the 
NACA 0015 airfoil at M=0.7. 



Figure 5-12 Effect of Pivot Point on reduced frequency for NACA 0015 and M=0.7 (Euler case) 


The area below the curve is an unstable region for the airfoil; that is if the 
frequency of oscillation is lower than the value corresponding to the curve for a certain 
pivot point, then flutter is triggered by an initial disturbance of 0.5 degrees. In contrast 
the area above the curve is the stable region and every oscillation will die out. 

It is readily seen that the most fluttei* prone range of pivot points is between 0 and 
-0.50-c. The pivot point most likely to induce flutter is Xp=-0.25-c, i.e when the pivot 
point of oscillation is one quarter chord ahead of the L.E. 

The above information can be shown in terms of the reduced natural pitching 
frequency ka required to avoid flutter. It should be noted that for elastic axis locations of 
-0.85<Xp<0.15 the used value of moment of inertia iawas 100. For Xp=0.15 and -0.85 


63 





the solutions were always decaying with this value of io so an increased value of 300 was 
used.The resulting graph is shown in Figure 5-13: 
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Figure 5-13 Effect of Pivot Point on critical kot for NACA 0015, M=0.7 and ia=100 (Euler case) 


Again, the area below the curve is an unstable region for the airfoil; if the airfoil 
spring constant Kais lower than the value that corresponds to the curve for a certain pivot 
point (low torsional stiffness), then flutter is induced by an initial disturbance of 0.5 rees. 

G. EFFECT OF MACH NUMBER ON TORSIONAL FLUTTER 

Figure 5-14 expands the information presented in Figure 5-12 by displaying the 
effect of Mach number on torsional flutter. 

An increase in Mach number significantly increases the region of instability for 
the whole range of pivot points. 

Another important result that comes from this graph is that the values of reduced 
frequency tend to become constant when reducing the Mach number to very low subsonic 
speeds of M=0.1 -0.2. 
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NACA 0015 



Xp 

Figure 5-14 Effect of Mach Number and Pivot Point for NACA 0015 (Euler case) 


This inforaiation is shown in bar chart form in Figure 5-15: 



Figure 5-15 Reduced frequency variation with Mach number for Pivot Points •0.85<Xp<0.15 
(NACA 0015) 
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The percentage increase of the reduced flutter frequency with Mach number and 
pivot points is tabulated in Table 5-4. 


M \ Xp 

0.15 

0 

-0.25 

-0.5 

-0.75 

-0.85 1 

0.7 

12.4% 

34.6% 

59.8% 

68.6% 

45.8% 


0.5 

11.4% 

13.8% 

23.4% 

28.9% 

28.3% 


0.3 

7.4% 

1.9% 

4.9% 

6.9% 

13.3% 


0.2 

2.4% 

-1.1% 

-2.1% 

-2.0% 

-0.1% 

4.0% 

0.1 

0 

0 

0 

0 

0 

0 


TaMe 5-4 Percentage Increase of k with Mach number and Pivot Points 


Another way to display the effect of pivot location and Mach number is to plot the 
reduced natural torsional frequency ka needed to suppress flutter. Again it is seen from 
Figure 5-16 that an increase in Mach number is destabilizing. 



Xp 


Figure 5-16 Effect of Mach number and pitch axis location on torsional flutter of NACA 0015 
airfoil, ia=100 
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Another useful way of representing the results is to show the variation of the 
reduced flutter velocity with Mach number for every pivot point. As discussed in Chapter 

IQ the reduced flutter velocity is defined as Vr = = — and can be used in order to 

0><xC K 

find the speed U for which flutter occurs for given values of Ka. Figure 5-17 shows the 
decrease of the flutter speed as the Mach number is increased. The results shown are only 
for pivot points -0.85<Xp<0 for which ia=100. 



Xp 


Figure 5-17 Variation of flatter speed with Mach Number for NACA 0015, i„=100 

H. EFFECT OF AIRFOIL THICKNESS 

The same procedure was applied to study the effect of airfoil thickness. 
Calculations were done for the airfoils: NACA 0006, NACA 0009, NACA 0012 and 
NACA 0015 for a Mach number of 0.7. The variation of critical reduced frequency with 
the airfoil thickness is presented in Figure 5-18 and Figure 5-19 showing that over the 
whole pivot point range any increase in airfoil thickness has a destabilizing effect. 
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Figure 5-18 Effect of airfoil thickness on critical reduced frequency at M=0.7 



Figure 5-19 Effect of airfoil thickness at M=0*7 
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This destabilizing influence of airfoil thickness is quantified in Table 5-5 where 
the percentage increase of the reduced flutter frequency k is tabulated. 


Airfoil 

thickness 

0.15 

0 

-0.25 

-0.5 

-0.75 

-0.85 








6 

0 

0 

0 

0 

0 

0 

9 

5.6% 

6.8% 

6.3% 



3.2% 

12 

11.2% 

37.8% 

49.8% 

18.8% 

5.2% 

5.4% 

15 

28.1% 

61.8% 

70.4% 

58.6% 

14.3% 

10.3% 


Table 5-5 Percentage increase of k with airfoU thickness and pivot point location 


Similarly, Figure 5-20 displays the reduced natural pitching frequency ka needed 
to suppress flutter. Again it is seen that the NACA 0015 airfoil requires a much higher 
torsional frequency than the NACA 0006 airfoil. Again the results shown are for elastic 



Xp 


Figure 5-20 Effect of airfoO thickness on torsional flutter at M=0.7, ia=100 
axis locations of -0.85<Xp<0.15 for which the used value of moment of inertia was 100. 
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Finally the destabilizing effect is also shown in Figure 5-21. Increasing the airfoil 
thickness results in a significant decrease in the flutter speed. 



Xp 


Figure 5-21 Variation of flutter speed with airfoil thickness at M=0.7, ia=100 


/. VISCOUS (NA VIER-STOKES) CASE CALCULATION 


The same procedure previously described for Euler case will be shown for the 
viscous flow over the same oscillating airfoil (NACA 0015), the same Mach number 
(M=0.7) and the same elastic axis of oscillation (L.E of the airfoil - pivot point Xp=0). 
The input file for the steady solution is as follows: 


# 

IREAD 

ITER 

NPRINT 

NLOAD 

ODVAR 



0 

20000 

1 

1 

O 

O 


# 

. ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 


0.5 

false 

false 

0.100 

-.500' 

0.500 

# 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 



false 

0. 

0.10 

0. 

0.1 

o 

o 

# 

MACH 

RE 

Vise 

TURBL 




0.700 

1.0e6 

true 

true 
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# 

TIMEAC 

COUR 

NEWTIT 






true 

1000.0 

1 





# 

free 

mass 

ialpha 

ka 

kh 

xp 

xa 


false 

0.0 

100.0 

0.00000 

0.0 

o 

o 

0.0 

# • 

aneut 

hneut 

hO 






0.0 

0.0 

0.0 






The steady solution is presented in Figure 5-22. We can see that Cl reaches a 
steady value of 0.08977 after t =40 sec. 
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Figure 5-22 Steady state soln for NACA 0015 M=0.7 (N-S case) 

The final result for the steady solution is shown at 
Table 5-6: 


T 

Aoa 

Cl 

Cd 

Cm 


0.5 

0.089777 

0.006215 

0.000665 


Table 5-6 Steady state solution for NACA 0015 M=0.7 (N-S case) 
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After the steady solution is obtained, the flutter condition can be found using the 
same procedure described for the Euler calculations. We start with a pivot point of 
Xp=0.0 and an assumed value of 0.36 for ka- The following is the input file for the N.S. 
code: 


# 

IREAD 

ITER 

NPRINT 

NLOAD 

ODVAR 



-1 

30000 

1 

1 

1.00 


# 

ALPHA 

OSCIL 

RAMP 

REDFRE 

ALFAMND 

ALFAMXD 


0.5 

false 

false 

0.100 

-.500 

0.500 

# 

PLUNGE 

PLMX 

PLMY 

PLPHSXD 

PLFREQ 



false 

0. 

0.10 

0. 

0.1 

0.0 0.0 

# 

MACH 

RE 

Vise 

TURBL 




0.700 

1.0e6 

true 

true 



# . 

TIMEAC 

COUR 

NEWTIT 





true 

700.00 

3 




# 

free 

mass 

ialpha 

ka 

kh 

xp xa 


true 

0.0 

100.0 

0.340000 

0.0 

0.0 0.0 

# 

aneut 

hneut 

hO 





0.0 

0.0 

0.0 





The plot of the resulting airfoil motion is presented in Figure 5-23. 



Figure 5-23 Unsteady case solution - Ka=0.34 (N-S case) 
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Since the airfoil motion is seen to decay the assumed value of ka is too high. After 
finding the maximum points of the graph and applying a regression equation we find that 
the slope of the maximum points is -3.188E-04 (Figure 5-24). 



Figure 5-24 Time rate of change of AOA amplitude for KopO.34 (N-S case) 
Similarly for ka=0.20 the diverging motion presented in Figure 5-25 is obtained. 













Then applying a regression equation for the maximum points of the graph, the 
slope is found to be 4.566E-04 (Figure 5-26): 


y = 4.383E-04X + 5.744E-01 



Figure 5-26 Time rate f hange of AOA amplitude for Ka=0.20 (N-S case) 


The results for the two trials are shown in Table 5-7 which make it possible to interpolate 
the value of ka which produces a constant amplitude of oscillation. As shown in Figure 5- 
27 this value is ka=0.252. 


Ka 

Slope 

0.34 

-2.860E-04 

0.20 

4.383E-04 


Table 5-7 Regression equation slope values for various ka values 
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y = -5.174236E-03X + 1.473201E-03 
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Figure 5-27 Curvefit of Ka vs time rate of change of AOA amplitude (N-S case) 


Using this ka value indeed produces the oscillation shown in Figure 5-28 and the 
slope -6.783E-10’^ shown in Figure 5-29. 
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Figure 5-28 Unsteady case solution - Ka=0.252 Critical flutter case (N-S) 
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y = -6.783E-05X + 5.466E-01 



time (NS) 


Figure 5-29 Time rate of change of AOA amplitude for Ka=0.324 Critical flutter (N-S case) 

From the airfoil motion we find that the average period of oscillation is 33.215 
(non-dimensional time) and reduced frequency of oscillation 0.270. 


J. EULER-NS RESULTS COMPARISON 

Due to time constraints the N-S study was performed only for NACA 0015 airfoil 
and M=0.7. The critical reduced frequency values for pivot points between 
-0.85<Xp<0.15 are presented in Figure 5-30 and compared with the Euler predictions. 

It is seen that the values of critical reduced frequencies found with the Euler 
calculations are up to 23% higher than the N.S predictions. As expected, viscous flow 
effects reduce the possibility of flutter. From Eq.(4.5) it is seen that the same airfoil at the 
same Mach number with 60% less torsional rigidity will have the same flutter stability. 
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Figure 5-30 Euler - NS reduced frequency results for NACA 0015 and M=0.7 


Table 5-8 shows the numerical differences between the above calculations for 
Euler and Navier-Stokes cases. 


Xp 

0.15 

0.00 

-0.25 

-0.50 

-0.75 

-0.85 








k(Euler) 

0.223 

0.339 

0.380 

0.334 

0.233 

0.221 

k (N-S) 

0.191 

0.270 

0.310 

0.260 

0.217 

0;205 

Difference 

14.3% 

20.3% 

18.4% 

22.1% 

6.8% 

7.2% 


Table 5-8 Numerical differences of reduced frequency between Euler and N-S calculations for 
NACA 0015 and M=0.7 


The reduced natural pitching frequency ka values computed with the N-S and 
Euler codes are presented in Figure 5-31 indicating that the minim u m torsional stiffness 
necessary to prevent flutter is smaller using viscous flow calculations. 
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In order to explain this finding it is instructive to examine the pressure distri¬ 
butions and Mach number contour plots for both cases in steady AOA of 0.5° (Figure 5- 
32 to Figure 5-39). The main difference of the graphs has to do with the strength of the 
shock above the airfoil. For the Euler case it is seen that a relatively strong shock is 
formed over the airfoil. The graphs that come from the Navier-Stokes calculations show 
that the shock over the airfoil is very weak. This is due to the boundary layer which 
smooths the shock. 



Xp 


Figure 5-31 Euler - NS reduced natural pitching frequency results for NACA 0015 and M=0.7 
4=100 
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Figure 5-32 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (Euler) 



Figure 5-33 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (Euler-detail) 
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Figure 5-34 Mach contours for NACA 0015 M=0.7 steady AOA (Euler) 



Figure 5-35 Mach contours for NACA 0015 M=0.7 steady AOA (Euler-detail) 







Figure 5-37 Pressure distribution contours for NACA 0015 M=0.7 steady AOA (N-S-detafl) 
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Figure 5-38 Mach contours for NACA 0015 M=0.7 steady AOA (N-S) 
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K. PANEL -EULER CODE COMPARISON 


Another interesting area of study is the comparison of results from the Euler code 
for non-viscous low subsonic flow with those from the UPOT panel code which 
calculates inviscid incompressible flow over the airfoil. The calculations were again 
made for flow the NACA 0015. The results for M=0.3,0.2,0.1 and incompressible flow 
are presented in Figure 5-40. 



NACA 0015 


Figure 5-40 Low subsonic Euler code results comparison with UPOT results. 

It is seen that the trends predicted by both codes are in good agreement, but there 
is a significant quantitative difference between the Euler prediction at M=0.1 and the 
incompressible panel code prediction. This difference could be reduced substantially if 
the panel code was run with a very small time step, as shown in Figure 5-40 and Figure 
5-41. 
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% difference 



Figure 5-41 % differences between Euler and UPOT results (w.r.t M=0.1 results) 
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VI. SUMMARY 


The time domain flutter analysis of NACA 0006, NACA 0009, NACA 0012 and 
NACA 0015 airfoils presented in this thesis leads to the following major conclusions: 

A. Linearized incompressible and linearized subsonic compressible flow 
theory yields unconservative estimates of single-degree-of-freedom torsional airfoil 
flutter. 

B. Torsional flutter region increases with increasing airfoil thickness and 
Mach number. 

C. Viscous flow effects have a stabilizing influence. 

D. Pitch axis locations upstream of the quarter chord point may induce flutter, 
especially in the range between the leading edge and a point a half-chord upstream of the 
leading edge. 

E. Flutter boundaries computed with the incompressible panel code UPOT 
require very small time steps; however, even for the smallest time step used, the UPOT 
computed flutter boundaries differed from Euler computed boundaries by a significant 
amount. 

F. Euler computations converged relatively quickly for high subsonic Mach 
numbers, but required increasingly large times for completion as the Mach number was 
•decreased toward 0.1. Every Euler computation could be completed in about 12 to 24 
hours. 

G. Navier-Stokes computation times typically were 4 to 6 times longer on the 
Department’s SGI work stations. 
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VII. RECOMMENDATIONS 


It will be interesting to study the effect of the following parameters: 

A. Reynolds Number 

Extend the calculation to Reynolds numbers other than 1-10^ 

B. Mach Number 

Extend the calculations to cover the transonic Mach number range and 
investigate the effect of shock motion and shock boundary layer interaction on 
flutter 

C. Moment of Inertia 

Investigate the effect of moment of inertia variation on flutter 

D. Two-Degree-of-Freedom Flutter 

Extend the single-degree-of-freedom analysis to the bending-torsion 
flutter problem and investigate the effect of airfoil thickness on this type of flutter. 


87 



88 



LIST OF REFERENCES 


1. Fung, Y.C, “Theory of Aeroelasticity,” Dover 1969. 

2. Fletcher, C.A., “Computational Techniques for Fluid Dynamics,” Springer 
Verlag,VOLI-n, 1990. 

3. Dowell, E.H., et al “A Modem Course in Aeroelasticity,” Sijthoff & Noordhoff 
1978. 

4. Bisplinghoff, R.L., and Ashley H., “Principles of Aeroelasticity,” Dover, 1962. 

5. Krothapalli A., “Recent Advances in Aerodynamics,” Springer Verlag, 1983. 

6. Bisplinghoff R.L., Ashley H., Halfman R.L., “Aeroelasticity,” Addison-Wesley, 
1951. 

7. Lamboume,N.C., “Flutter in One Degree of Freedom” Part V Chapter 5, AGARD 
Manual on Aeroelasticity, General Ed. W. P. Jones, Oct 1961. 

8. Schlichting,H., “Boundary Layer Theory, McGraw-Hill Book Company,”, 1966. 

9. Cricelli, A.M., “Unsteady Airfoil Flow Solutions on Moving Zonal Grids,” 
Master’s Thesis, Naval Postgraduate School, 1992. 

10. Baldwin, B.S., Lomax, H., “Thin Layer Approximation and Algebraic Model for 
Separated Turbulent Flows,” AIAA Paper 78-257, Jan 1978. 

11. Rai, M.M. and Chakravarthy, S.R., “An Implicit Form of the Osher Upwind 
Scheme,” AIAA Journal Vol.24, No 5, May 1986, pp-735-743. 

12. Cowles L.J., “High Reynolds Number, Low Mach Number Steady Flow Field 
Calculations over NACA 0012 Airfoil Using Navier-Stokes and Interactive 
Boundary Layer Theory,” Master’s Thesis, Naval Postgraduate School, 1987. 

13. Jones K.D. and Platzer M.F., “Airfoil Geometry and Flow Compressibility Effects 
on Wing and Blade Flutter,” AIAA Paper 98-0517, January 1998. 

14. Jones K.D. and Platzer M.F., “Time Domain Analysis of Low-Speed Airfoil 
Flutter,” AIAA Journal, Vol 34, No 5, May 1996. 


89 



15. Teng, N.H., “The Development of a Computer Code for the Numerical Solution 
of Unsteady, Inviscid and Incompressible Flow over an Airfoil,” Master’s Thesis, 
Naval Postgraduate School, 1987. 


90 



INITIAL DISTRIBUTION LIST 


1. Defense Technical Information Center.2 

8725 John J.Kingman Road, Ste 0944 

Ft.Belvoir, VA 22060-6218 

2. Dudley Knox Library .2 

Naval Postgraduate School 

411 Dyer Rd. 

Monterey, CA 93943-5101 

3. Chairman.1 

Department of Aeronautics and Astronautics, Code AA 

Naval Postgraduate School 
699 Dyer Rd, Room 137 
Monterey, CA 93943-5106 

4. Dr. Max F.Platzer.7 

Department of Aeronautics and Astronautics, Code AA/PL 

Naval Postgraduate School 
699 Dyer Rd, Room 137 
Monterey, CA 93943-5106 

5. Dr. Kevin. D. Jones.1 

Department of Aeronautics and Astronautics, Code AA 

Naval Postgraduate School 
699 Dyer Rd, Room 137 
Monterey, CA 93943-5106 


91 








6. CAPT. Constantines Kakkavas 
Peloponnisou 19 
Vrilissia 15235 
Athens, GREECE 




